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Several methods for preparing well equilibrated melts of long chains polymers are studied. We 
show that the standard method in which one starts with an ensemble of chains with the correct 
end-to-end distance arranged randomly in the simulation cell and introduces the excluded volume 
rapidly, leads to deformation on short length scales. This deformation is strongest for long chains 
and relaxes only after the chains have moved their own size. Two methods are shown to overcome 
this local deformation of the chains. One method is to first pre-pack the Gaussian chains, which 
reduces the density fluctuations in the system, followed by a gradual introduction of the excluded 
volume. The second method is a double-pivot algorithm in which new bonds are formed across a pair 
of chains, creating two new chains each substantially different from the original. We demonstrate 
the effectiveness of these methods for a linear bead spring polymer model with both zero and 
nonzero bending stiffness, however the methods are applicable to more complex architectures such 
as branched and star polymer. 
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I. INTRODUCTION 

A systematic investigation of the structure-property 
relations for polymeric systems by computer simulations 
requires the preparation of equilibrated melts of long, 
entangled chains. For temperatures well above the glass 
transition, this can, in principle, be achieved using suf- 
ficiently long molecular dynamics (MD) or Monte Carlo 
(MC) simulations. However since the longest relaxation 
of an entangled polymer melt of length N scales at least 
as N 3 , giving at least A 4 in cpu time, this method is 
only feasible for relatively short chain lengths. Depend- 
ing on the detail of the model the longest chains which 
can presently be equilibrated by brute-force MD or MC 
simulations are on the order of 2-7 entanglements lengths 
N e . Using united atom potentials for polyethylene in 
which one treats the carbon and its bonded hydrogen as a 
single united atom, it is possible to equilibrate high tem- 
perature (T > 45QK) melts for chains of the order of ap- 
proximately 2N e or about 200 monomers. 1,2 For coarse- 
grained bead-spring models like the one employed here 
it is possible to study longer chain lengths, on the order 
of up to 500 monomers or approximately 7N e . 3 However 
this reaches the very limits of present day fastest comput- 
ers. An increase of cpu speed by an order of magnitude 
does even not allow for a doubling of the chain length. 
While the present chain lengths are sufficient to follow 
the dynamics well into the reptation regime, they are 
not long enough to study structure-property relations, 
such as the plateau modulus G° N . This would require 
chains of order 10 — 20A e . The situation is even worse 
for branched polymer melts or polymers at interfaces, 
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where equilibration times even for relatively short chains 
are prohibitively large to use brute force simulations to 
produce equilibrated melts. 

Fortunately to produce equilibrated samples, there is 
no need to follow the slow physical dynamics of entan- 
gled polymer melts. One possibility is to temporarily 
turn off the excluded volume interactions and to allow 
the chains to pass through each other. In order to ob- 
tain a well-defined Monte Carlo algorithm, it is useful to 
combine this idea with parallel tempering. While it has 
been demonstrated that this is feasible in principle, the 
method has so far not proven particularly efficient, the 
main reason being the large amount of computer time 
spent on unphysical Hamiltonians. More conventional 
MC algorithms which can be used to equilibrate poly- 
mer melts include reptation moves (generalized slithering 
snake algorithms), 5 configuration bias algorithms, 6-8 and 
concerted rotation algorithms. 8-10 Being relatively lo- 
cal, these methods work best for moderate chain lengths 
and densities. The complete equilibration of very long 
chain melts still requires long runs. The fastest method, 
the generalized slithering snake algorithm, theoretically 
scales as A 2 . An alternative are algorithms, which are 
able to move large sections of a chain at once. The 
prototype for such methods is the pivot algorithm 11-15 
in which a monomer is chosen at random and one of 
the two segments formed by that partitioning is piv- 
oted rigidly in a random direction about the unit. A 
Boltzmann weight is used to determine if the move is ac- 
cepted or not. This method is highly efficient for single 
chains in a continuum, implicit solvent. Unfortunately, 
a direct application to dense melts is not feasible since 
any large scale conformational change is bound to vio- 
late the packing constraints imposed by the chain ex- 
cluded volume. However it is possible to introduce a 
double bond crossing algorithm in a way that two new 
bonds or bridges arc formed across a pair of chains, creat- 
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ing two new chains each substantially different from the 
original. 2 ' 10 ' 16-18 For the bead spring model concerned 
here this can simply be done by cutting two bonds and 
introducing two new bonds as discussed in more detail 
in Sec. VII. For atomistic models, the bond length is 
much shorter then the diameter of a monomer, and it 
is necessary to construct new bridges between the chain 
segments. Karayiannis et al. 17,ls have shown that by 
constructing two trimer bridges between the two prop- 
erly chosen dimers along the backbone one can quickly 
equilibrate long chain polyethylene melts. For the two 
chains involved this is essentially a double-pivot move, 
where each chain experiences such a pivot rotation. 

In addition to improved equilibration algorithms there 
is an obvious interest in methods for generating ini- 
tial melt configurations which are as close to equilib- 
rium as possible. 19-24 In our previous work on polymer 
melts, 3 ' 19 ' 25,26 we prepared most melts by first generat- 
ing an ensemble of chains with the correct end-to-end 
distance and randomly placing them in the simulation 
cell. We then introduced a weak, non-diverging excluded 
volume potential in the form of e.g. a cosine potential, 
A(l + cos(7rr/2 1 / 6 er)), between non-bonded monomers, 
where r is the distance between two monomers and a is 
the bead diameter. The amplitude A was then increased 
over a short time interval until the minimum distance be- 
tween monomers was sufficient to switch on the Lennard- 
Jones potential without creating numerical instabilities. 
In the first part of this paper we confirm observations 
by Brown et al. 21 that this method deforms the (longer) 
chains. As a consequence, long chain melts are not as 
well equilibrated as originally believed. The deformation 
turns out to be strongest at short length scales and com- 
pletely relaxes away only after a time T m ax(N), which is 
the time for a chain to move its own size. 21 This long 
time is needed because the local chain-chain topology 
can only relax by the slow physical dynamics. For lin- 
ear chains, T max {N) « r e (N/N e ) 3 , where the entangle- 
ment time r e = TR. OUS e(iV e ) is the Rouse time of a chain 
of length N e (experimentally the measured exponent is 
closer to 3.4 than 3 except for extremely long chains). 
This however was exactly what we tried to avoid by that 
approach. Not surprisingly, the longer the chain length 
N, the more severe the effect is. For relatively short 
chains, such as those investigated in our previous stud- 
ies of the crossover from Rouse to rcptation dynamics 
(N < 350), the simulations were run long enough that the 
results are independent of the starting state. 3 ' 19 ' 25 How- 
ever for longer chains, which cannot be run long enough 
compared to the longest relaxation time, this simple pro- 
cedure for generating starting states is inadequate. 

In the present paper we report results from an ex- 
tensive effort to prepare well-equilibrated melts of bead- 
spring polymers with chain lengths ranging from N = 
350 to N — 7000. For this purpose we use both ap- 
proaches outlined above. First we show how to avoid 
the local stretching with only minimal computational ef- 
fort by a suitable modification of our standard method. 



Then we demonstrate the capacity of the double-pivot 
algorithm to dynamically equilibrate melts of medium 
sized chains of length up to N = 700. The paper is or- 
ganized as follows: In Sec. II we define the model. In 
Sec. Ill, we present some simple theoretical estimates of 
the structure of bead-spring chains with different intrinsic 
bending stiffness in dense melts. In Sec. IV, we discuss 
ways to characterize the single chain structure and ex- 
tract suitable target functions for long chain melts from 
long simulations of short chain melts. In Sec. V we ex- 
amine carefully the standard procedure used in the past 
to prepare melt configurations. We show that for long 
chains this method deforms the chains at short length 
scales and that these deformations relax completely only 
after the chains have moved their own size. In Sec. VI, 
we present the new pre-packing procedure, which signif- 
icantly reduces the density fluctuations particularly for 
long chains. This combined with a gradual introduc- 
tion of the excluded volume, results in well equilibrated 
long chain melts in a reasonable amount of cpu time. In 
Sec. VI we describe the double-pivot algorithm and com- 
pare results for the internal dimensions of the chain with 
those produced from pre-packing and a gradual introduc- 
tion of the excluded volume. In Sec. VIII we test the two 
methods for chains with local bending rigidity. A brief 
summary of our main conclusions are given in Sec. IX. 



II. MODEL 

We use a coarse grained model in which the polymer 
is treated as a string of beads of mass m connected by 
a spring. The beads interact with a purely repulsive 
Lcnnard-Jones excluded volume interaction, 
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in addition to the Lennard-Joncs interaction. The model 
parameters are the same as in ref. 19, namely k = 30e/(T 2 
and R a = 1.5a unless otherwise noted. The temperature 
T = e/kB- The basic unit of time is r = a(mj 'e) 1 / 2 . We 
performed constant volume simulations of monodisperse 
melts at a segment density p = 0.85er~ 3 . The temper- 
ature was kept constant by coupling the motion of each 
bead weakly to a heat bath with a local friction T. Un- 
less otherwise specified T = 0.5m/r. The equations of 
motion were integrated using a velocity Vcrlet algorithm 
with a time step At = 0.012r. The average bond length 
is < b 2 > 1 / 2 = 0.97cr. The polymer melts studied con- 
sisted of M chains of length TV beads. The chain lengths 
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studied varied from N — 25 to 7000. The number of 
chains are each system is specified in Table I. 

Equilibration algorithms for polymer melts which fol- 
low the physical dynamics encounter a strong growth of 
the necessary equilibration times with chain length. Per 
unit of time r, the simulation of a melt of M chains of 
length N requires a computational effort of 10 2 A/ x N 
particle updates. Since the longest relaxation time is of 
the order r max rj r e (N/N e ) 3 , equilibration of our largest 
system with M = 80 and N = 7000 would require 
10 2 M x N{T e /r)(N/N e ) 3 w 8 x 10 17 particle updates. 
Assuming a typical performance of 2.5 x 10 5 particle up- 
dates per second, which is typical for this model on a 
single processor, the required cpu time for the equilibra- 
tion is about 10 5 years. On a modern parallel cluster such 
as the DEC alpha CPlant cluster at Sandia, the times are 
somewhat less. On 64 processors, we get about 20 steps 
per second for our largest system, which means that the 
total time for the chains to move on average their own 
size is approximately 1600 years. Increasing the numbers 
of processors helps somewhat but since the efficiency of 
the computation decreases as the number of processors 
is increased for fixed system size, this type of brute force 
approach for long chain melts is clearly not feasible. 

In our previous studies, we considered only fully flex- 
ible chains. However to demonstrate the effectiveness of 
the methodology, we also include a nonzero bending stiff- 
ness, which is modelled by 
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where cos 6^ = (?< - • (f m - f<), where (f, - r<_i) 
is the unit vector in the bond direction. 



III. CHAIN CHARACTERISTIC IN THE MELT 

In dense polymer melts excluded volume interactions 
between different parts of a given polymer chain are 
screened. Since in our model there are no explicit in- 
trachain interactions beyond those between neighboring 
bonds, the single chain structure should be well charac- 
terized by the expectation value (cos 9) of the bond angle 
8. For example, the knowledge of (cos 9) is sufficient to 
calculate an estimate for the mean-square end-to-end ex- 
tension (R 2 ) of chain segments of length N monomers 
and n = N — 1 bonds, 



{R 2 )(n)=nb 2 



'l + (cos<9) 12(cos6>)(l-(cos0) n )' 



1- (cos 6) n (1- (cos 6>})' 

where the asymptotic prefactor is referred to as 

1 + (cos 9) 



1 — (cos ( 



FIG. 1. Estimates of the effective stiffness Coo of our chains 
as a function of the strength kg of the bending potential 
Eq. (2). Eq. (5) (solid line) only accounts for the bending 
energy, while Eq. (6) (•) also considers excluded volume in- 
teractions between next-nearest neighbor monomers along the 
chain. The gray dots indicate estimates based on the analy- 
sis of simulation data for equilibrated melts of chain length 
N = 200 - 350 for kg = and 50 - 100 for kg > 0. 



The simplest estimate of (cos(6*)) only takes the bare 
bending energy Eq. (2) into account. Neglecting small 
variations in the bond length around the mean value b = 
0.97cr, one finds: 
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where (3 — 1/ksT. Eq. (5) is represented as a solid 
line in Fig. 1. However this result underestimates the 
effective chain stiffness, since the chain cannot fold back. 
This effect can be accounted for approximately by mod- 
elling the chains as non-reversal-random- walks (NRRWs) 
with excluded volume interactions between next-nearest 
neighbors along the backbone. As a function of the bond 
angle 9 their distance is given by lg = — r"i-i\ = 

6[2(l + cos6» i )] 1 / 2 so that 
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which has to be evaluated numerically (see the black dots 
in Fig. 1). For kg = this gives Coo = 1.76 compared to 
Cqo = 1 obtained from Eq.( III) and in good agreement 
with simulation results for chains of length 20 < N < 
350. This value is slightly smaller than the value c^, = 
2.06 obtained at the O-point for the same Lennard- Jones 
interaction truncated at 2.5<r in an implicit, continuum 
solvent. 15 

Ensembles of single chains which obey Eq. (3) can be 
generated by simple sampling. This is particularly sim- 
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pie, if the allowed bond vectors are distributed evenly 
over a solid angle with 6 < max . In this case 
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so that the appropriate value for 9 max can be easily de- 
termined for any desired effective stiffness. 



IV. TARGET FUNCTIONS AND SUITABLE 
WAYS OF CHARACTERIZING SINGLE CHAIN 
STRUCTURE 

The estimates of the chain structure presented in the 
previous section were obtained within an effective sin- 
gle chain picture. In reality, the formal characterization 
of a polymer melt is a complicated many-body problem. 
For example, we have completely ignored the influence of 
packing effects on the chain conformations. As was re- 
cently shown 27 the local melt structure very sensitively 
depends on the ratio of bond length to effective excluded 
volume of the beads. This effect is not well represented by 
the value of the overall chain extension. Here we use our 
standard bond length/diameter ratio and focus on the 
single chain properties, however the methods discussed 
below apply to systems with different parameters as well. 
In the present case, sufficiently long MD or MC simula- 
tions represent the only "ab initio" method which allows 
to take into account all interactions. 

In order to characterize the chain conformations in 
the melt we mainly rely on mean-square internal dis- 
tances (R 2 ) N) averaged over all segments of size 
n = \i — j\ along the chains, where i < j G [1, -W] are 
monomer indices. In addition we characterize deviations 
from Gaussian statistics on short length scales using the 
ratio (i? 4 ) / (R 2 ) ■ Figures 1 and 2 show a compari- 
son of simulation results for equilibrated melts and our 
naive estimates. The results from the melts have been 
run sufficiently long that the chains have moved several 
times their own size. For N = 350, kg — 0, the run time 
T t = 2.6 x 10 6 r. 3 Clearly, the description of chains in a 
melt as simple freely jointed chains is too naive. Nev- 
ertheless, the deviations are not large and the measured 
and estimated values for Coo agree quite well. In particu- 
lar, the simulation data show neither unexpected features 
nor significant finite chain length effects. 

The rest of the present paper deals with algorithms 
which try to circumvent the slow physical equilibration 
path. These algorithms will be judged according to two 
criteria: their capacity to reproduce the target functions 
and their performance. We mostly concentrate on the 
fully flexible case. In the end, we come back to the stiffcr 
systems as a kind of "blind test" . 
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FIG. 2. Mean square internal distances from long MD runs 
for chain length N = 100 and 350 respectively for four values 
of the intrinsic bending stiffness: kg = (□), kg — 0.75 (x), 
kg = 1.5 (*), and kg = 2.0 (+). We refer to these data sets as 
"target functions" . For comparison, they are included as thick 
black lines in subsequent, corresponding plots. The dashed 
lines show to a fit of Eq. (3) to the asymptotic behavior. The 
corresponding values for are included in Fig. 1. 



V. STANDARD PROCEDURE FOR PREPARING 
POLYMER MELTS 

Our standard method for preparing melt conforma- 
tions is as follows: 

1. Start from an ensemble of chains with the correct 
end-to-end distance R 2 (N) = 6 2 c 00 (A r — 1) on large 
length scales. For a known Coo this step is trivial. 

2. Arrange the chains randomly in the simulation box. 

3. Introduce a weak, non-diverging excluded volume 
potential. A convenient form for the soft potential 
is a cosine potential, 
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between non-connected monomers with 
The initial amplitude of A = 4e was linearly in- 
creased to a final values of A > lOOe ("push-off") 
over a short time interval of 10-20 r. In the final 
state the inter monomer distances are large enough 
to allow one to switch to the LJ potential without 
creating numerical instabilities. During this phase, 
we often use a stronger coupling to the thermostat 
by increasing the friction constant to T = 2.0to/t. 
In addition, the velocities of all monomers are set to 
zero every 50 time steps. This facilitates monomers 
that are strongly overlapping to separate. We re- 
fer to this method in which the excluded volume 
interactions are introduced fairly rapidly as "fast 
push-off." 
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FIG. 3. Mean square internal distances after a fast push-off 
for randomly packed phantom chains of length N — 25 (+), 
50 (x), 350 (*) and 7000(D). The thick line is our target 
function from Fig. 2. 



FIG. 4. Mean square internal distances for the systems 
shown in Fig. 3 after an additional MD relation of 10 4 r which 
corresponds to tr ousc (N = 90). Symbols as in Fig. 3. The 
thick line is our target function from Fig. 2. 



4. Relax the system with the full LJ-potential by a 
short MD run. 

Figure 3 shows that this procedure actually deforms 
the chains so that the melt is not equilibrated after step 
3. The deformations are strongest on short lengths. 
With increasing TV only the amplitude but not the po- 
sition of the deformation along the chain changes. The 
monomer displacements during step 3 are too small to 
affect the conformations on large scales, allowing us to 
tune them to preselected values. Since the largest length 
scales also take the longest time to relax, one might hope 
that the MD relaxation in step 4 can be significantly 
shorter than for a starting conformation of, say, com- 
pletely stretched chains. To be more specific, a naive 
expectation is that chain segments of length n < N equi- 
librate on time scales comparable to the Rouse/reptation 
time of chains of length n independent of N. Fig- 
ure 4 compares the result of equilibration runs of length 
T t = 10 4 r « T max (N = 90) to the target function in 
Fig. 2. As expected, the internal distances measured for 
chains of length N = 25 and N = 50 coincide. However 
segments with n < 50 of longer chains for N = 350 arc 
still far from being equilibrated after 10 4 r. Even after 
a time t ps tr ousc (N) ~ T e (N/N e ) 2 ~ 10 5 r (not shown) 
the chains of length N = 350 were not fully equilibrated. 
Rather the local equilibration is completed only after the 
chains have moved their own size, t > T max (N). In lat- 
tice simulations it is sometimes possible to avoid the slow 
reptation dynamics by allowing the chains to cut through 
each other. 28,29 In this case T max is given by a Rouse- 
time and scales as N 2 . In off-lattice bead-spring models 
topology conservation is the result of energetic barriers 
which prevent chain crossing. Since these barriers are a 
result of the microscopic interaction potentials, they are 
difficult to manipulate without affecting the local chain 
structure. Hence Rouse-like relaxation mechanisms due 
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500 


100 


250 


200 


120,200 


350 


500 


500 


200 


700 


200 


1400 


80 


3500 


80 


7000 



TABLE I. System size studied: Number of chains M and 
chain length N. 



to barrier crossing dominate reptation only for extremely 
long (and therefore inaccessible) chain lengths. Attempts 
to circumvent these barriers using parallel tempering 4 
have met with limited success at least as far as efficiency 
gains are concerned. 

VI. OPTIMIZED PROCEDURE FOR POLYMER 
MELT CONFORMATIONS 

The methods presented in this section aim at a more 
careful implementation of the idea underlying our stan- 
dard procedure: the local build-up of the characteristic 
melt packing for chains with the correct large length scale 
statistics. In this section we show that we can achieve 
this goal by first pre-packing the phantom chains and a 
subsequent slow and improved push-off scheme for the 
excluded volume. Both steps are needed to produce well 
equilibrated melts in a reasonable amount of cpu time. 
Applying only one does not achieve this objective. All 
results in this section are averaged over five independent 
realizations of the packing and push-off procedures. The 
system sizes studied are listed in Table I. 
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A. Prepacking of phantom chains 
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The chain deformations created by the standard pro- 
cedure are due in part to the large local density fluctua- 
tions in a "melt" of randomly overlapping NRRW chains. 
As the chain length increases, these fluctuations increase 
in size as shown in Fig. 5. In this section we describe 
a dynamic Monte-Carlo-like pre-packing procedure for 
phantom chains, which significantly reduces these den- 
sity fluctuations. 

As a merit function we have chosen the fluctuations 
E = (n 2 ) — (n) 2 in the number m of particles found in 
a sphere of radius d around particle i. In a homoge- 
neously dense system E = 0. We perform a zero temper- 
ature Monte Carlo (MC) optimization in which all moves 
which decrease the density fluctuations are accepted and 
all moves which increase the density fluctuations are re- 
jected. In the course of the packing run d is decreased 
from d = 4cr to 2a. We use a linked cell structure in order 
to calculate the rii efficiently. Choosing d > a effectively 
increases the range of the excluded volume correlations. 

All MC moves change the positions or orientations of 
entire chains which are treated as rigid. This pre-packing 
procedure therefore does not affect the single chain statis- 
tics, which by construction already have the correct end- 
to-end distance R as well as the targeted internal distance 
distributions. We consider five types of MC moves: 

Translation of individual chains in a random directions. 

Rotation of individual chains by random angles around 
random axes through their centers of mass. 

Reflection of individual chains at random mirror planes 
going through the center of mass. 

Inversion of individual chains at their centers of mass. 

Exchange of two chains preserving the center of mass 
positions. 

Typical run times were of the order of days on individual 
workstations and therefore negligible. 

The amplitude of the density fluctuations is drastically 
reduced as the packing proceeds. The final result is best 
assessed by the q — > limit of the system structure func- 
tion S(q). Figure 5 shows that the pre-packing reduces 
the density fluctuations to one percent of the initial value. 
Applying the fast push-off procedure outlined in Sec. V to 
these pre-packed states, we see by comparing Figs. 3 and 
6a that the local deformations arc eliminated for short 
chains (N < 50) but not for long chains. For long chains, 
there is a significant reduction in the local stretching but 
not enough to avoid the necessity of subsequent long MD 
runs, which again would need to be of the order of the 
disentanglement time T max (N). 
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FIG. 5. Density fluctuation of randomly distributed (□) 
and pre-packed Gaussian chains (■) for various chain lengths. 
By pre-packing the density fluctuations are about 1% for 
N = 350 and 0.5% for N = 7000 of that for randomly dis- 
tributed chains. 



B. Slow push-off 

In order to further reduce the chain deformations we 
have modified the push-off procedure in three ways: 

1. Replace the cos-potential Eq. (3) with a force- 
capped-Lennard- Jones-potential, 21 . 

2. Keep the full excluded volume between next- 
nearest neighbor monomers to preserve to non- 
reversal-random- walk structure (see sec. III). 

3. Increase the push-off time. 

Force-capping 21 can either be defined directly via a 
maximum force F max or implicitly via a distance rf c 
where Ulj (r/ c ) = F max . At larger separations, the in- 
teraction is defined by Eq. (1), at smaller distances the 
force becomes constant so that the interaction potential 
grows only linearly with (r — ff c ): 
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In the present case, rj c is gradually reduced from the 
Lennard-Jones cut-off radius r c = 2 x l % a to 0.8a which is 
significantly smaller than the relevant intcrparticle dis- 
tances. 

Force-capping has the advantage that the soft potential 
systematically approaches the true potential. In contrast, 
the cos-potential overcompensates the missing singular- 
ity at the origin by a fast rise of the repulsive interactions 
close to the cut-off distance. This leads to slight differ- 
ences in the monomer packing. As a result, the chain 
conformations are slightly perturbed after switching to 
the full LJ-potential. 
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push off procedure cannot be "too slow" and the chain 
size varies only very weakly with the maximum excluded 
volume force. 

Figure 6b shows that the new procedure works quite 
well. Independent of the total chain length N the chains 
arc no longer stretched on intermediate scales and the 
large scales remain unaffected by the introduction of ex- 
cluded volume interactions. As an additional check, we 
have varied the stiffness of the initially generated random 
walks over a range of about 20% (Fig. 7a). Fig. 7b com- 
pares the internal distances after pre-packing and slow 
push-off and confirms our previous observations. The 
chains arc fully equilibrated on short scales, where the 
results are independent of the initial conditions. In con- 
trast, internal distances beyond n = 100 remain prac- 
tically unaffected. We therefore conclude that our new 
method allows the preparation of well-equilibrated melts 
of extremely long chains within reasonable computational 
effort, provided the correct asymptotic chain stiffness Cqq 
is known from a careful extrapolation of simulation re- 
sults for well-equilibrated short and medium chain length 
melts. For our largest system of M = 80 chains of length 
N = 7000, an equilibration time of 5000r takes about 
one month on a single processor or about 8 hours on 
the 100 processor DEC alpha cluster. Thus even with 
moderate computational means it is possible to equili- 
brate very large systems using a combined pre-packing 
and slow push-off procedure. 



FIG. 6. Mean square internal distances for pre-packed sys- 
tems after a (a) fast push-off and (b) slow push-off for chains 
with kg = 0. Symbols as in Fig. 3. The thick line is our 
target function from Fig. 2. 



We have increased the push-off time in order to intro- 
duce the excluded volume interactions quasi-statically. 
In practice, one has in this context to worry about two 
issues: (i) How slow is slow enough? (ii) What is the 
equilibrium statistics for polymers with force-capped in- 
teractions? Concerning the first point, we have chosen a 
push-off time of 5000r, which is of the order tr O usc(50). 
This covers the typical subchain regime, up to which the 
intrachain distances still were disturbed by a fast push- 
off. The second question is more difficult to answer. If 
we simply switch off the repulsive monomer-monomer in- 
teraction, we find (-R 2 ) = 0.94cr 2 A r , while the effective 
stiffness Cqo ~ 1.7 of our chains in the melt is due to 
local bead packing effects. Thus a slow push off start- 
ing from a potential for which the end-to-end distance 
is an ideal random walk is also dangerous. For bead- 
spring chains these effects can easily be accounted for 
by treating the chains as NRRWs by always keeping the 
full LJ excluded volume between next-nearest neighbors 
along the chain. For atomistic polymer models, it is nec- 
essary to extend this to larger distances along the chain 
("pentane effect"). 21 Once this adjustement is made the 



VII. DOUBLE-PIVOT-MD HYBRID 
ALGORITHM 

The present section deals with the problem of acce- 
larating the dynamic equilibration of a long chain poly- 
mer melt. We describe a double-pivot-MD hybrid algo- 
rithm which performs this task significantly faster than 
simple MD relaxation. In contrast to the methods dis- 
cussed in the previous section, the dynamic equilibration 
of a dense polymer melt does not require independent 
knowledge of the effective Cqo. However, the preparation 
of locally equilibrated samples can significantly reduce 
the required relaxation times. 

In the original pivot algorithm 11-14 one choses a 
monomer at random and one of the two segments formed 
by that partitioning is pivoted rigidly in a random direc- 
tion about the selected monomer. The pivot algorithm 
is extremely efficient when applied to isolated chains in 
an implicit solvent. However, in a melt or even in a solu- 
tion at moderate density, pivot moves applied to a single 
chain are bound to be rejected, since they lead to strong 
excluded volume interactions with other chains. 

An attractive alternative are Monte Carlo moves in- 
volving at least two chains, which change the connectivity 
within the melt in a way that the overall packing of beads 
remains (almost) unaffected. Such algorithms were first 
introduced for lattice models 30,31 and recently extended 
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FIG. 7. Influence of the effective stiffness used to generate 
the initial sets of random walks on the final configurations for 
systems with ke = 0. The figure shows mean square inter- 
nal distances (a) before and (b) after pre-packing and slow 
push-off. From top to bottom the curves represent results for 
chains with length N = 3500 which were set up with a stiffness 
varying by -4%, -2%, 0%, +2%, +4%, +6%, +8%, +10% 
and +12% relative to our initial guess of Coo = 1.7. The thick 
line is our target function from Fig. 2. 

to off-lattice models. 2 ' 10 ' 16-18 The efficient equilibration 
usually comes at the expense of a certain amount of poly- 
dispcrsity which is introduced into the samples, though 
recently Karayiannis 17 ' 18 have overcome this limitation 
using double-bridging Monte Carlo moves. For the bead- 
spring polymers studied here, the algorithm is straight- 
forward to implement since the bond length I is compara- 
ble to the excluded volume parameter a. This certainly 
is a unique situation, as for many bead spring coarse 
grained models of specific chemical species, the above ra- 
tio is not as close to one. For example in the case of 
a united atom model for polyethylene where the bond 
length I = 1.54 A is significantly smaller then a w 4.0A, 
it is necessary to construct two trimer bridges between 
the two properly chosen dimers along the backbone. 17 ' 18 
For the bead spring model studied here, complex bridg- 
ing moves are not necessary. For this reason we refer 



to the method as the double-pivot algorithm since it is 
much closer to the pivot algorithm used for single chains 
in an implict, continuum solvent. 

In order to maintain a monodisperse melt we search 
for a pair of spatial neighbor monomers i and j on differ- 
ent chains or on the same chain, which happen to also be 
the same distance from one end of their respective chain. 
Then one can replace two bonds along the original chains 
with two new bonds or bridges across the pair of chains. 
The total change in energy AE is local and consists of 
the sum of the energies of the two new bonds minus the 
energies of the two previous bonds as well as the dif- 
ference in the bending energy, which involves the sum 
of the energies of four new angles minus the energies of 
the four previous angles for the case kg ^ 0. The move 
is accepted by a standard Metropolis criterion, namely 
with a Boltzmann weight exp (— AE/UbT) if AE > 
and 1 if AE < 0. As in the original pivot algorithm, the 
double-pivot algorithm generates new chains which are 
substantially different from the original two chains since 
as many as N/2 monomers of a chain are replaced by 
monomers from the other chain. This results in immedi- 
ate large changes in the end-to-end distance and radius 
of gyration. 

We have implemented the double-pivot algorithm into 
both our shared memory and distributed memory MD 
codes. The details of the implementation and a discus- 
sion of its computational efficiency will be published else- 
where. Briefly in our shared memory code, every 3 — 5 
time steps, we randomly chose a monomer i and check 
its non-bonded nearest neighbors to determine if any 
satisfy the condition that they are equidistant from the 
end of their chain. If so, the energy change AE of the 
double-pivot move is determined and the move accepted 
or rejected on the basis of the Metropolis criterion. If 
the move is accepted, monomers and/or the connectiv- 
ity table are re-labeled depending on which of the two 
codes is used and the MD simulation is continued. If the 
move is rejected or no suitable pair is identified, a new 
monomer is chosen at random and the process repeated. 
If no suitable pair is generated after searching a specified 
fraction of monomers (usually 2 — 5% depending on the 
chain length N), we continue with the MD simulations 
for another 3 — 5 steps and repeat the procedure. On 
the distributed memory code, the procedure is similar 
except that each processor randomly choses a monomer 
i and checks its non-bonded nearest neighbors on the 
same processor to see if they satisfy the condition that 
they are equidistant from the end of their chain. If no 
suitable non-bonded neighbor is found, another monomer 
i is randomly selected. The process is continued until a 
specified fraction of the momomers, typically 50%, are 
tested. To facilitate determining the distance from the 
free end, each chain carries an extra label starting from 
1 at either end to N/2 at the center. Because each pro- 
cessor searches a unique region of space independently, 
it is possible to have multiple successful moves each time 
the procedure is applied. The search is restricted to sets 
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of monomers on the same processor to avoid the need to 
communicate changes in chain topologies between proces- 
sors. This restriction also prevents two or more proces- 
sors from performing simultaneous swaps that could en- 
ergetically conflict with each other. While this restriction 
means (slightly) less swaps are considered, this is more 
than outweighed by the parallelism, e.g. up to P swaps 
take place in one time step, where P is the number of 
processors. With the standard parameters of the FENE 
potential, k = 30e/a 2 , the acceptance rate is quite low. 
A way to improve the acceptance rate is to start with a 
somewhat lower value, k = 10 — 15e/cr 2 and gradually in- 
crease k during the course of the simulation. This change 
has little effect on the overall structure of the chains but 
increases significantly the acceptance rate. Using our dis- 
tributed memory code on 27 DEC alpha processors for a 
system of M = 500 chains of length N = 500, the num- 
ber of successful exchanges for 10 5 Ai with kg = was 
approximately 1750 for k = 30e/<7 2 compared to 8500 for 
k = 20e/a 2 and 49500 for k = 10e/cr 2 . 

To demonstrate the effectiveness of the double-pivot al- 
gorithm, we have applied it to a melt of M = 500 chains 
of length N = 500 generated using the standard fast 
push-off procedure outlined in Sec. V. Figure 8 shows 
the internal distances as a function of time during a 
doublc-pivot/MD relaxation. Note that the character- 
istic "bump" in the internal distance function does not 
simply decay. Rather, Figure 8 shows that the (seem- 
ingly equilibrated) largest scales are first swollen, before 
the chain dimensions decay homogeneously on all scales. 
Equilibration for the reduced spring constant k = 10e/cr 2 
was achieved after about 6 x 10 4 r or 10 accepted pivot 
moves per monomer. The penalty for reducing the spring 
constant to k = 10 is an increase of the mean-square bond 
length by 7%, while the chain stiffness Coo remains un- 
affected. In subsequent runs, the spring constant k is 
increased to slowly "reel in" in the chains. This required 
an additional 6 x 10 4 r. We made no attempt to optimize 
the number of steps for each value of k. Nevertheless, 
the required equilibration times arc considerably shorter 
than the time r„ la2 -(500) ~ 4 x 10 6 r for the chains to 
move their own size in a pure MD simulation. 

These observations illustrate the close analogy between 
the pivot algorithm for single chains and the double-pivot 
algorithm for dense polymer melts. As pointed out by 
Madras and Sokal 13 , it is useful to distinguish between 
decorrelation and equilibration in discussing the perfor- 
mance of such algorithms. Pivoting is extremely efficient 
in decorrclating large length scales, provided the initial 
conformation is properly equilibrated. Up to small cor- 
rections, the number of pivot moves needed to decorrelate 
the end-to-end distance of isolated chains in an implicit 
solvent is O(10) independent of chain length. Since the 
same holds for subchains of arbitrary length, correspond- 
ingly more moves are required to decorrelate higher or- 
der (Rouse) modes. In other words, the decorrelation 
proceeds from large to small scales. However, the same 
is not correct for the equilibration of a perturbed ini- 



tial configuration. In this case it is necessary to run the 
system up to the decorrelation time of the shortest per- 
turbed length scale in order to equilibrate the chains on 
all length scales. For single chain studies aiming at global 
properties Madras and Sokal 13 therefore recommend to 
apply the pivot algorithm to equilibrated initial config- 
urations generated by other techniques. In the present 
case, the combination of pre-packing with a slow push- 
off at least partially fulfills these requirements, since the 
chains are locally equilibrated. 

Since equilibration requires that each monomer comes 
close enough to suitable exchange partners, one may ask 
whether of O(10) exchanges per monomer independent of 
N is sufficient or whether the exchanges are simply being 
made back and forth between a limited number of states. 
To estimate whether enough independent configurations 
are being sampled, consider the time for a monomer 
to diffusive a typical distance between monomers with 
the same index (p/(N/2))~ 1 / 3 between successful pivot 
moves. If we assume that monomers diffuse with the 
typical (5r 2 ) = ^/kBTb 2 t/( Rouse dynamics, this corre- 
sponds to 



tpiv 



k R Tb 2 



(8) 



which is equivalent to the Rouse time of a chain of length 
N ef f = {N/2) 2 / 3 . For our longest chains N = 7000, 
N e ff w 236. The Rouse time for a chain of length 
TV = 236 is approximately 7 x 10 4 t, which is still much 
less than the time needed to make O(10) exchanges per 
monomer even if we used the reduced spring constant 
k = 10. Thus only for much longer chains or more com- 
plex architectures does one have to worry about sufficient 
independent configurations being sampled by this hybrid 
doublc-pivot-MD algorithm. 



VIII. LOCAL BENDING RIGIDITY 

As a last point we present results for chains with in- 
trinsic stiffness Eq. (2). In section IV we referred to this 
case as a "blind test," because the starting states were 
generated based on our simple estimate Eq. (6) for the 
effective chain stiffness in the melt and before the brute- 
force MD runs for the target functions shown in Fig. 2 
were completed. As can be seen in Fig. 1, our original 
estimates were slightly too large. A comparison with the 
target functions after the pre-packing and slow push-off 
phase (Fig. 9a) shows that as for fully flexible chains 
the correct chain statistics is reproduced on short length 
scales. However, we are now in a situation comparable to 
the situtaion presented in Fig. 7 where the large length 
scales are not fully equilibrated, because the initial chains 
were generated with a slightly incorrect effective stiffness. 

Subsequently, we ran the combined double-pivot-MD 
simulation for 1.08 x 10 5 t with k = 10 for the first 4 
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FIG. 8. Equilibration of mean square internal distances 
using combined double-pivot-MD relaxation following a fast 
push-off for 500 chains of length N = 500. Results for 
t = (+) just after the push-off, 6 x 10 3 r (x), 1.2 x 10 3 r (*), 
6 x 10 4 r (□) and 1.2 x 10 5 r (■). In order to increase the 
acceptance rate of the pivot moves we used a reduced spring 
constant fc = 10e/a 2 for < t < 6 x 10 4 r and k = 20e/a 2 for 
6 x 10 4 r < t < 9.6 x 10 4 r before switching to fc = 30e/V 2 for 
9.6 x 10 4 r < t < 1.2 X 10 5 t. 





fc = 10 


fc = 20 


fc = 30 


kg = 0.0 


24800 


3900 


840 


kg = 0.75 


27000 


4200 


894 


kg = 1.5 


21200 


3500 


770 


kg = 2.0 


15500 


2700 


600 



TABLE II. Influence of the strength of the FENE spring 
constant fc on the number of successful double-pivot exchanges 
per 10 5 time steps for a system of 200 chains of length 
N = 350 on 27 DEC alpha processors. 



millions steps, fc = 20 for the second 4 million steps and 
fc = 30 for the last million steps for a system of 200 
chains of N = 350 for kg = 1.5 and 2.0 and N = 700 for 
kg = 0.75. The resulting mean squared end-to-end dis- 
tance were in excellent agreement with the target func- 
tions generated by brute force MD simulations for shorter 
chains. Table II shows an example of the number of 
successful exchanges per 10 5 steps for the system of 200 
chains of length N = 350 monomers for three values of 
the spring constant k for four values of kg for our dis- 
tributed memory code run on 27 DEC alpha processors. 
As seen from this table, the procedure used gives approx- 
imately 9 — 12 successful exchanges for k = 10, 2 success- 
ful exhanges for fc = 10 and 0.1 for fc = 30. Fortunately 
there is little difference in the local structure of the chain 
for fc = 20 compared to fc = 30, though as seen from 
Fig. 8, the chain is expanded for fc = 10. Thus there is a 
definite tradeoff between reducing the bond spring con- 
stant fc, which increases the number of exchanges while at 
the same time increasing the mean squared bond length. 
Note that for longer chains lengths, the length of the 
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FIG. 9. Mean square internal distances for chains with lo- 
cal bending rigidity (a) after pre-packing and slow push-off 
and (b) after an additional double-pivot-MD relaxation of 
length 1.08 x 10 5 r. Results for kg = (*), 0.75e (+), 
1.5e (x), and 2.0e (□). We did not apply the additional dou- 
ble-pivot-MD relaxation to our fully flexible systems, which 
are well-equilibrated after pre-packing and slow push-off. The 
thick lines are our target functions from Fig. 2 which were 
generated a posteriori for chains with intrinsic stiffness. 

run and/or the number of processors would have to be 
increased so that the number of successful exchanges is 
on the order of O(10). This limits the double-pivot/MD 
hydrid method to moderate chain lengths. 



IX. CONCLUSIONS 

In this paper, we have discussed the preparation of 
equilibrated melts of long chain polymers in computer 
simulations. The interest in reliable algorithms for this 
task is due to two problems: (i) the prohibitively long 
relaxation times encountered in brute force MD simula- 
tions and (ii) the difficulty to predict the chain statistics 
on large scales (i.e. the effective chain stiffness Coo) on the 
basis of microscopic intra- and interchain interactions. 

Our results can be summarized as follows: 
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• It is insufficient to judge the quality of the equi- 
libration from the statistics of the chain end-to- 
end distances or radius of gyration. The first just 
measures one length, while the latter gives a non 
specific average over all internal distances. As a 
much more sensitive criterion we measure mean 
square internal distances on all length scales from 
the monomer size up to the contour length of the 
chains. 

• Our standard method to rapidly introduce ex- 
cluded volume interactions between randomly as- 
sembled Gaussian chains with the correct overall 
statistics ("fast push-off") fails, when judged ac- 
cording to this improved criterion. The fast push- 
off introduces significant chain distortions on length 
scales of the order of the tube diameter. They can 
only vanish when the local melt topology is equili- 
brated. In MD simulations this is only possible via 
reptation dynamics so that the proper equilibration 
requires runs whose length exceeds the disentangle- 
ment time of the chains. To overcome this problem 
we have introduced two different methods. 

• We have modified our standard approach by re- 
ducing the density fluctuations in the assembly of 
Gaussian chains ( "pre-packing" ) and by introduc- 
ing the excluded volume interactions in a quasi- 
static manner ("slow push-off"). As before, this 
method requires prior knowledge of Coo- Tests 
for bead-spring models with chain lengths up to 
N = 7000 (N/N e = 0(100)) show that suitable 
push-off times are the Rouse times of the character- 
istic chain lengths of the overshoot. In the present 
case this time is of the order of a few entanglement 
times, independent of N. 

• We have applied a "double-pivot" algorithm which 
is inspired by the highly efficient pivot algorithm 
for single chains and the double-bridging method 
for dense systems. For intermediate chain lengths, 
this is a viable method to equilibrate melts in a 
reasonable amount of cpu time, particularly when 
one has no prior knowledge of Cqo. 

The combination of MD relaxation, double-pivot and 
slow push-off allows the efficient and controlled prepa- 
ration of equilibrated melts of short, medium and long 
chains respectively. While our results were obtained for 
an off-lattice bead spring model with chain lengths up 
to N = 7000 beads, the methods should work as ef- 
ficiently for lattice polymer models such as the bond 
fluctuation model. 32,33 The pre-packing and gradual in- 
troduction of the excluded volume are also applicable 
to united and explicit atom simulations. Furthermore, 
the methods are applicable to polydisperse melts and 
to branched polymers, including long chain branching, 
comb and star polymers. The double-pivot and double 
bridging algorithms has also been used to equilibrate long 



end-tethered, polymer brushes in contact with a melt of 
long chains. 18,34 For these more complex chain architec- 
tures, in which the chains are not necessarily Gaussian 
and there is no a priori way to know the conformation of 
the chains, algorithms like the double-pivot and double- 
bridging algorithms are presently the only possible way 
to equilibrate systems faster than extremely long MD or 
MC runs. 
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